##############################################################
## Replication code
## Clarke, "Revolutionary Violence and Counterrevolution"
##############################################################

## Packages ###########################################################################################################
library(brglm)
library(lmtest)
library(sandwich)
library(margins)
library(ggplot2)
library(extrafont)

## Data ###########################################################################################################
## Load data 
CR <- read.csv("counterrev.csv")

## Create data subset of revs that witness counterrev challenges
CRsub <- CR[CR$CRbin == 1,]

## Main models ###########################################################################################################
## Deaths
deaths <- brglm(sucCRbin ~ lndeaths 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                + incumbpowerdur + incgovmilitary
                + leftist + vanguard + nr_spons
                ,
                data = CR, family = binomial)

deaths_cl <- coeftest(deaths, vcov = vcovCL, cluster = ~ cowcode)

deathschal <- brglm(CRbin ~ lndeaths 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbpowerdur + incgovmilitary
                    + leftist + vanguard + nr_spons
                    ,
                    data = CR, family = binomial)

deathschal_cl <- coeftest(deathschal, vcov = vcovCL, cluster = ~ cowcode)

deathssuc <- brglm(sucCRbin ~  lndeaths 
                   + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                   + incumbpowerdur + incgovmilitary
                   + leftist + vanguard + nr_spons
                   
                   ,
                   data = CRsub, family = binomial)

deathssuc_cl <- coeftest(deathssuc, vcov = vcovCL, cluster = ~ cowcode)

## Civil war
cw <- brglm(sucCRbin ~ civilwar 
            + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
            + incumbpowerdur + incgovmilitary
            + leftist + vanguard + nr_spons
            ,
            data = CR, family = binomial)

cw_cl <- coeftest(cw, vcov = vcovCL, cluster = ~ cowcode)

cwchal <- brglm(CRbin ~ civilwar 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                + incumbpowerdur + incgovmilitary
                + leftist + vanguard + nr_spons
                ,
                data = CR, family = binomial)

cwchal_cl <- coeftest(cwchal, vcov = vcovCL, cluster = ~ cowcode)

cwsuc <- brglm(sucCRbin ~  civilwar 
               + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
               + incumbpowerdur + incgovmilitary
               + leftist + vanguard + nr_spons
               ,
               data = CRsub, family = binomial)

cwsuc_cl <- coeftest(cwsuc, vcov = vcovCL, cluster = ~ cowcode)

## Rev militia
revmil <- brglm(sucCRbin ~ revcoerc 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                + incumbpowerdur + incgovmilitary
                + leftist + vanguard + nr_spons
                ,
                data = CR, family = binomial)

revmil_cl <- coeftest(revmil, vcov = vcovCL, cluster = ~ cowcode)

revmilchal <- brglm(CRbin ~ revcoerc 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbpowerdur + incgovmilitary
                    + leftist + vanguard + nr_spons
                    ,
                    data = CR, family = binomial)

revmilchal_cl <- coeftest(revmilchal, vcov = vcovCL, cluster = ~ cowcode)

revmilsuc <- brglm(sucCRbin ~  revcoerc 
                   + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                   + incumbpowerdur + incgovmilitary
                   + leftist + vanguard + nr_spons
                   ,
                   data = CRsub, family = binomial)

revmilsuc_cl <- coeftest(revmilsuc, vcov = vcovCL, cluster = ~ cowcode)

## Figures #################################################################################################
## Deaths figures
plotdata_deaths <- cplot(deaths, "lndeaths", vcov = vcovCL(deaths, cluster = ~ cowcode), draw = FALSE)
plotdata_deathschal <- cplot(deathschal, "lndeaths", vcov = vcovCL(deathschal, cluster = ~ cowcode), draw = FALSE)
plotdata_deathssuc <- cplot(deathssuc, "lndeaths", vcov = vcovCL(deathssuc, cluster = ~ cowcode), draw = FALSE)

plot_deaths <- ggplot() +
  geom_line(data = plotdata_deaths, aes(x = xvals, y = yvals)) +
  geom_ribbon(data = plotdata_deaths, aes(x = xvals, ymin = lower, ymax = upper), alpha = .2) +
  theme_bw(base_family = "Garamond") +
  ylim(c(-0.1,1)) +
  ylab("Probability of Counterrevolution") +
  xlab("Deaths (log)") +
  ggtitle("Counterrevolution (aggregate outcome)") +
  theme(axis.text.x = element_text(size = 14, color = "black")) +
  theme(axis.text.y = element_text(size = 14, color = "black")) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_text(size = 14, color = "black")) +
  theme(plot.title = element_text(hjust = 0.5, size = 16, color = "black"))
ggsave("plot_deaths.png", plot = plot_deaths, width = 5, height = 4.5)

plot_deathschal <- ggplot() +
  geom_line(data = plotdata_deathschal, aes(x = xvals, y = yvals)) +
  geom_ribbon(data = plotdata_deathschal, aes(x = xvals, ymin = lower, ymax = upper), alpha = .2) +
  theme_bw(base_family = "Garamond") +
  ylim(c(-0.1,1)) +
  ylab("Probability of Counterrev challenge") +
  xlab("Deaths (log)") +
  ggtitle("Counterrevolutionary challenge") +
  theme(axis.text.x = element_text(size = 14, color = "black")) +
  theme(axis.text.y = element_text(size = 14, color = "black")) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_text(size = 14, color = "black")) +
  theme(plot.title = element_text(hjust = 0.5, size = 16, color = "black"))
ggsave("plot_deathschal.png", plot = plot_deathschal, width = 5, height = 4.5)

plot_deathssuc <- ggplot() +
  geom_line(data = plotdata_deathssuc, aes(x = xvals, y = yvals)) +
  geom_ribbon(data = plotdata_deathssuc, aes(x = xvals, ymin = lower, ymax = upper), alpha = .2) +
  theme_bw(base_family = "Garamond") +
  ylim(c(-0.1,1)) +
  ylab("Probability of Counterrev success") +
  xlab("Deaths (log)") +
  ggtitle("Counterrevolutionary success") +
  theme(axis.text.x = element_text(size = 14, color = "black")) +
  theme(axis.text.y = element_text(size = 14, color = "black")) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_text(size = 14, color = "black")) +
  theme(plot.title = element_text(hjust = 0.5, size = 16, color = "black"))
ggsave("plot_deathssuc.png", plot = plot_deathssuc, width = 5, height = 4.5)

## Binary vars figures
vc <- vcovCL(revmil, cluster = ~ cowcode)
colnames(vc) <- labels(coefficients(revmil))
rownames(vc) <- labels(coefficients(revmil))
marg_revmil <- margins(revmil, variable = c("revcoerc"), vcov = vc)
marg_revmil_t <- summary(marg_revmil)

vc <- vcovCL(cw, cluster = ~ cowcode)
colnames(vc) <- labels(coefficients(cw))
rownames(vc) <- labels(coefficients(cw))
marg_cw <- margins(cw, variable = c("civilwar"), vcov = vc)
marg_cw_t <- summary(marg_cw)

plotdata_binary <- matrix(NA, ncol = 3, nrow = 2)
colnames(plotdata_binary) <- c("prob", "lower", "upper") 
rownames(plotdata_binary) <- c("Rev Militia", "Civil War")

plotdata_binary[1,1] <- marg_revmil_t[,2]
plotdata_binary[1,2] <- marg_revmil_t[,6]
plotdata_binary[1,3] <- marg_revmil_t[,7]

plotdata_binary[2,1] <- marg_cw_t[,2]
plotdata_binary[2,2] <- marg_cw_t[,6]
plotdata_binary[2,3] <- marg_cw_t[,7]

plotdata_binary <- as.data.frame(plotdata_binary)
plot <- ggplot(plotdata_binary, aes(colour = rownames(plotdata_binary)))
plot_binary <- plot + geom_linerange(aes(x = factor(rownames(plotdata_binary), 
                                                         levels = c("Rev Militia", "Civil War")), 
                                              y = prob, ymin = lower, ymax = upper), 
                                          size = 1.5, position = position_dodge(width = 4/5)) +
  geom_pointrange(aes(x = factor(rownames(plotdata_binary),
                                 levels = c("Rev Militia", "Civil War")),
                      y = prob, ymin = lower, ymax = upper),
                  fatten = 4, linetype = "solid", size = 1, shape = 20) +
  ylab("Marginal Effect") + xlab("") + theme_bw(base_family = "Garamond") + 
  theme(legend.position = "none") + ylim(-0.5,0.5) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14, color = "black")) +
  theme(axis.text.y = element_text(size = 14, color = "black")) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_text(size = 14, color = "black")) + coord_flip() + 
  ggtitle("Counterrevolution (aggregate outcome)") +
  theme(plot.title = element_text(hjust = 0.5, size = 16, color = "black")) +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, color = "dark grey") +
  scale_color_manual(name = "", values = c("darkolivegreen3", "tomato3"))

png("plot_binary.png", width = 4, height = 4, units = 'in', res = 300)
plot_binary
dev.off()

vc <- vcovCL(revmilchal, cluster = ~ cowcode)
colnames(vc) <- labels(coefficients(revmilchal))
rownames(vc) <- labels(coefficients(revmilchal))
marg_revmilchal <- margins(revmilchal, variable = c("revcoerc"), vcov = vc)
marg_revmilchal_t <- summary(marg_revmilchal)

vc <- vcovCL(cwchal, cluster = ~ cowcode)
colnames(vc) <- labels(coefficients(cwchal))
rownames(vc) <- labels(coefficients(cwchal))
marg_cwchal <- margins(cwchal, variable = c("civilwar"), vcov = vc)
marg_cwchal_t <- summary(marg_cwchal)

plotdata_binarychal <- matrix(NA, ncol = 3, nrow = 2)
colnames(plotdata_binarychal) <- c("prob", "lower", "upper") 
rownames(plotdata_binarychal) <- c("Rev Militia", "Civil War")

plotdata_binarychal[1,1] <- marg_revmilchal_t[,2]
plotdata_binarychal[1,2] <- marg_revmilchal_t[,6]
plotdata_binarychal[1,3] <- marg_revmilchal_t[,7]

plotdata_binarychal[2,1] <- marg_cwchal_t[,2]
plotdata_binarychal[2,2] <- marg_cwchal_t[,6]
plotdata_binarychal[2,3] <- marg_cwchal_t[,7]

plotdata_binarychal <- as.data.frame(plotdata_binarychal)
plot <- ggplot(plotdata_binarychal, aes(colour = rownames(plotdata_binarychal)))
plot_binarychal <- plot + geom_linerange(aes(x = factor(rownames(plotdata_binarychal), 
                                                    levels = c("Rev Militia", "Civil War")), 
                                         y = prob, ymin = lower, ymax = upper), 
                                     size = 1.5, position = position_dodge(width = 4/5)) +
  geom_pointrange(aes(x = factor(rownames(plotdata_binarychal),
                                 levels = c("Rev Militia", "Civil War")),
                      y = prob, ymin = lower, ymax = upper),
                  fatten = 4, linetype = "solid", size = 1, shape = 20) +
  ylab("Marginal Effect") + xlab("") + theme_bw(base_family = "Garamond") + 
  theme(legend.position = "none") + ylim(-0.5,0.5) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14, color = "black")) +
  theme(axis.text.y = element_blank()) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_blank()) + coord_flip() + 
  ggtitle("Counterrevolutionary challenge") +
  theme(plot.title = element_text(hjust = 0.5, size = 16, color = "black")) +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, color = "dark grey") +
  scale_color_manual(name = "", values = c("darkolivegreen3", "tomato3"))

png("plot_binarychal.png", width = 3, height = 4, units = 'in', res = 300)
plot_binarychal
dev.off()

vc <- vcovCL(revmilsuc, cluster = ~ cowcode)
colnames(vc) <- labels(coefficients(revmilsuc))
rownames(vc) <- labels(coefficients(revmilsuc))
marg_revmilsuc <- margins(revmilsuc, variable = c("revcoerc"), vcov = vc)
marg_revmilsuc_t <- summary(marg_revmilsuc)

vc <- vcovCL(cwsuc, cluster = ~ cowcode)
colnames(vc) <- labels(coefficients(cwsuc))
rownames(vc) <- labels(coefficients(cwsuc))
marg_cwsuc <- margins(cwsuc, variable = c("civilwar"), vcov = vc)
marg_cwsuc_t <- summary(marg_cwsuc)

plotdata_binarysuc <- matrix(NA, ncol = 3, nrow = 2)
colnames(plotdata_binarysuc) <- c("prob", "lower", "upper") 
rownames(plotdata_binarysuc) <- c("Rev Militia", "Civil War")

plotdata_binarysuc[1,1] <- marg_revmilsuc_t[,2]
plotdata_binarysuc[1,2] <- marg_revmilsuc_t[,6]
plotdata_binarysuc[1,3] <- marg_revmilsuc_t[,7]

plotdata_binarysuc[2,1] <- marg_cwsuc_t[,2]
plotdata_binarysuc[2,2] <- marg_cwsuc_t[,6]
plotdata_binarysuc[2,3] <- marg_cwsuc_t[,7]

plotdata_binarysuc <- as.data.frame(plotdata_binarysuc)
plot <- ggplot(plotdata_binarysuc, aes(colour = rownames(plotdata_binarysuc)))
plot_binarysuc <- plot + geom_linerange(aes(x = factor(rownames(plotdata_binarysuc), 
                                                        levels = c("Rev Militia", "Civil War")), 
                                             y = prob, ymin = lower, ymax = upper), 
                                         size = 1.5, position = position_dodge(width = 4/5)) +
  geom_pointrange(aes(x = factor(rownames(plotdata_binarysuc),
                                 levels = c("Rev Militia", "Civil War")),
                      y = prob, ymin = lower, ymax = upper),
                  fatten = 4, linetype = "solid", size = 1, shape = 20) +
  ylab("Marginal Effect") + xlab("") + theme_bw(base_family = "Garamond") + 
  theme(legend.position = "none") + ylim(-0.5,0.5) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 14, color = "black")) +
  theme(axis.text.y = element_blank()) +
  theme(axis.title.x = element_text(size = 14, color = "black")) +
  theme(axis.title.y = element_blank()) + coord_flip() + 
  ggtitle("Counterrevolutionary success") +
  theme(plot.title = element_text(hjust = 0.5, size = 16, color = "black")) +
  geom_hline(yintercept = 0, linetype = "dashed", size = 1, color = "dark grey") +
  scale_color_manual(name = "", values = c("darkolivegreen3", "tomato3"))

png("plot_binarysuc.png", width = 3, height = 4, units = 'in', res = 300)
plot_binarysuc
dev.off()

## Descriptive statistics #################################################################################################
sumstats <- matrix(NA, ncol = 5, nrow = 15)
colnames(sumstats) <- c("Mean", "SD", "Min", "Max", "Source")
rownames(sumstats) <-  c("Counterrevolution",
                         "Number of deaths in revolution (log)", "Civil war", "Revolutionary militia",
                         "GDP per capita (log)", "Population (log)","Urbanization %", "Ethnic fractionalization",
                         "% territory mountainous (log)","End year of revolution",
                         "Duration incumbent was in office","Incumbent military regime",
                         "Leftist ideology in revolution", "Revolution headed by vanguard party", 
                         "Rev regime client of major world power")
sumstats[,5] <- c("Original data", 
                  "Beissinger 2022", "Beissinger 2022", "Original data",
                  "Maddison","EUGene","EUGene and Vanhannen","Wimmer and Min 2006","Fearon and Laitin 2003",
                  "Calculated","Beissinger 2022","Beissinger 2022",
                  "Beissinger 2022", "Beissinger 2022", "Casey 2020")

sumstats[1,1] <- mean(na.omit(CR$sucCRbin))
sumstats[1,2] <- sd(na.omit(CR$sucCRbin))
sumstats[1,3] <- min(na.omit(CR$sucCRbin))
sumstats[1,4] <- max(na.omit(CR$sucCRbin))

sumstats[2,1] <- mean(na.omit(CR$lndeaths))
sumstats[2,2] <- sd(na.omit(CR$lndeaths))
sumstats[2,3] <- min(na.omit(CR$lndeaths))
sumstats[2,4] <- max(na.omit(CR$lndeaths))

sumstats[3,1] <- mean(na.omit(CR$civilwar))
sumstats[3,2] <- sd(na.omit(CR$civilwar))
sumstats[3,3] <- min(na.omit(CR$civilwar))
sumstats[3,4] <- max(na.omit(CR$civilwar))

sumstats[4,1] <- mean(na.omit(CR$revcoerc))
sumstats[4,2] <- sd(na.omit(CR$revcoerc))
sumstats[4,3] <- min(na.omit(CR$revcoerc))
sumstats[4,4] <- max(na.omit(CR$revcoerc))

sumstats[5,1] <- mean(na.omit(CR$lngdpcapplus1))
sumstats[5,2] <- sd(na.omit(CR$lngdpcapplus1))
sumstats[5,3] <- min(na.omit(CR$lngdpcapplus1))
sumstats[5,4] <- max(na.omit(CR$lngdpcapplus1))

sumstats[6,1] <- mean(na.omit(CR$lntpop))
sumstats[6,2] <- sd(na.omit(CR$lntpop))
sumstats[6,3] <- min(na.omit(CR$lntpop))
sumstats[6,4] <- max(na.omit(CR$lntpop))

sumstats[7,1] <- mean(na.omit(CR$urbpercaftrev))
sumstats[7,2] <- sd(na.omit(CR$urbpercaftrev))
sumstats[7,3] <- min(na.omit(CR$urbpercaftrev))
sumstats[7,4] <- max(na.omit(CR$urbpercaftrev))

sumstats[8,1] <- mean(na.omit(CR$ethfrac))
sumstats[8,2] <- sd(na.omit(CR$ethfrac))
sumstats[8,3] <- min(na.omit(CR$ethfrac))
sumstats[8,4] <- max(na.omit(CR$ethfrac))

sumstats[9,1] <- mean(na.omit(CR$lnmtnest))
sumstats[9,2] <- sd(na.omit(CR$lnmtnest))
sumstats[9,3] <- min(na.omit(CR$lnmtnest))
sumstats[9,4] <- max(na.omit(CR$lnmtnest))

sumstats[10,1] <- mean(na.omit(CR$endyear))
sumstats[10,2] <- sd(na.omit(CR$endyear))
sumstats[10,3] <- min(na.omit(CR$endyear))
sumstats[10,4] <- max(na.omit(CR$endyear))

sumstats[11,1] <- mean(na.omit(CR$incumbpowerdur))
sumstats[11,2] <- sd(na.omit(CR$incumbpowerdur))
sumstats[11,3] <- min(na.omit(CR$incumbpowerdur))
sumstats[11,4] <- max(na.omit(CR$incumbpowerdur))

sumstats[12,1] <- mean(na.omit(CR$incgovmilitary))
sumstats[12,2] <- sd(na.omit(CR$incgovmilitary))
sumstats[12,3] <- min(na.omit(CR$incgovmilitary))
sumstats[12,4] <- max(na.omit(CR$incgovmilitary))

sumstats[13,1] <- mean(na.omit(CR$leftist))
sumstats[13,2] <- sd(na.omit(CR$leftist))
sumstats[13,3] <- min(na.omit(CR$leftist))
sumstats[13,4] <- max(na.omit(CR$leftist))

sumstats[14,1] <- mean(na.omit(CR$vanguard))
sumstats[14,2] <- sd(na.omit(CR$vanguard))
sumstats[14,3] <- min(na.omit(CR$vanguard))
sumstats[14,4] <- max(na.omit(CR$vanguard))

sumstats[15,1] <- mean(na.omit(CR$nr_spons))
sumstats[15,2] <- sd(na.omit(CR$nr_spons))
sumstats[15,3] <- min(na.omit(CR$nr_spons))
sumstats[15,4] <- max(na.omit(CR$nr_spons))

sumstats[,1] <- round(as.numeric(sumstats[,1]), digits = 2)
sumstats[,2] <- round(as.numeric(sumstats[,2]), digits = 2) 
sumstats[,3] <- round(as.numeric(sumstats[,3]), digits = 2)
sumstats[,4] <- round(as.numeric(sumstats[,4]), digits = 2)

## Counterrevolutionary counts ######################################################################################
deathschal_CRcount <- glm(CRchallenge ~ lndeaths 
                          + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                          + incumbpowerdur + incgovmilitary
                          + leftist + vanguard + nr_spons
                          ,
                          data = CR, family = poisson)

deathschal_CRcount_cl <- coeftest(deathschal_CRcount, vcov = vcovCL, cluster = ~ cowcode)

cwchal_CRcount <- glm(CRchallenge ~ civilwar 
                      + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                      + incumbpowerdur + incgovmilitary
                      + leftist + vanguard + nr_spons
                      ,
                      data = CR, family = poisson)

cwchal_CRcount_cl <- coeftest(cwchal_CRcount, vcov = vcovCL, cluster = ~ cowcode)

revmilchal_CRcount <- glm(CRchallenge ~ revcoerc 
                          + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                          + incumbpowerdur + incgovmilitary
                          + leftist + vanguard + nr_spons
                          ,
                          data = CR, family = poisson)

revmilchal_CRcount_cl <- coeftest(revmilchal_CRcount, vcov = vcovCL, cluster = ~ cowcode)

deathssuc_CRcount <- brglm(sucCRbin ~  lndeaths + CRchallenge
                           + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                           + incumbpowerdur + incgovmilitary
                           + leftist + vanguard + nr_spons
                           
                           ,
                           data = CRsub, family = binomial)

deathssuc_CRcount_cl <- coeftest(deathssuc_CRcount, vcov = vcovCL, cluster = ~ cowcode)

cwsuc_CRcount <- brglm(sucCRbin ~  civilwar + CRchallenge
                       + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                       + incumbpowerdur + incgovmilitary
                       + leftist + vanguard + nr_spons
                       ,
                       data = CRsub, family = binomial)

cwsuc_CRcount_cl <- coeftest(cwsuc_CRcount, vcov = vcovCL, cluster = ~ cowcode)

revmilsuc_CRcount <- brglm(sucCRbin ~  revcoerc + CRchallenge
                           + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                           + incumbpowerdur + incgovmilitary
                           + leftist + vanguard + nr_spons
                           ,
                           data = CRsub, family = binomial)

revmilsuc_CRcount_cl <- coeftest(revmilsuc_CRcount, vcov = vcovCL, cluster = ~ cowcode)

## Robustness test #############################################################################################
## Add controls incrementally ####################################################################################
## Deaths
deaths1 <- brglm(sucCRbin ~ lndeaths
                 ,
                 data = CR, family = binomial)

deaths_cl1 <- coeftest(deaths1, vcov = vcovCL, cluster = ~ cowcode)

deaths2 <- brglm(sucCRbin ~ lndeaths 
                 + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                 ,
                 data = CR, family = binomial)

deaths_cl2 <- coeftest(deaths2, vcov = vcovCL, cluster = ~ cowcode)

deaths3 <- brglm(sucCRbin ~ lndeaths 
                 + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                 + incumbpowerdur + incgovmilitary
                 ,
                 data = CR, family = binomial)

deaths_cl3 <- coeftest(deaths3, vcov = vcovCL, cluster = ~ cowcode)

deathschal1 <- brglm(CRbin ~ lndeaths 
                     ,
                     data = CR, family = binomial)

deathschal_cl1 <- coeftest(deathschal1, vcov = vcovCL, cluster = ~ cowcode)

deathschal2 <- brglm(CRbin ~ lndeaths 
                     + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                     ,
                     data = CR, family = binomial)

deathschal_cl2 <- coeftest(deathschal2, vcov = vcovCL, cluster = ~ cowcode)

deathschal3 <- brglm(CRbin ~ lndeaths 
                     + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                     + incumbpowerdur + incgovmilitary
                     ,
                     data = CR, family = binomial)

deathschal_cl3 <- coeftest(deathschal3, vcov = vcovCL, cluster = ~ cowcode)

deathssuc1 <- brglm(sucCRbin ~  lndeaths 
                    ,
                    data = CRsub, family = binomial)

deathssuc_cl1 <- coeftest(deathssuc1, vcov = vcovCL, cluster = ~ cowcode)

deathssuc2 <- brglm(sucCRbin ~  lndeaths 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    ,
                    data = CRsub, family = binomial)

deathssuc_cl2 <- coeftest(deathssuc2, vcov = vcovCL, cluster = ~ cowcode)

deathssuc3 <- brglm(sucCRbin ~  lndeaths 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbpowerdur + incgovmilitary
                    ,
                    data = CRsub, family = binomial)

deathssuc_cl3 <- coeftest(deathssuc3, vcov = vcovCL, cluster = ~ cowcode)

## Civil war
cw1 <- brglm(sucCRbin ~ civilwar 
             ,
             data = CR, family = binomial)

cw_cl1 <- coeftest(cw1, vcov = vcovCL, cluster = ~ cowcode)

cw2 <- brglm(sucCRbin ~ civilwar 
             + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
             ,
             data = CR, family = binomial)

cw_cl2 <- coeftest(cw2, vcov = vcovCL, cluster = ~ cowcode)

cw3 <- brglm(sucCRbin ~ civilwar 
             + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
             + incumbpowerdur + incgovmilitary
             ,
             data = CR, family = binomial)

cw_cl3 <- coeftest(cw3, vcov = vcovCL, cluster = ~ cowcode)

cwchal1 <- brglm(CRbin ~ civilwar 
                 ,
                 data = CR, family = binomial)

cwchal_cl1 <- coeftest(cwchal1, vcov = vcovCL, cluster = ~ cowcode)

cwchal2 <- brglm(CRbin ~ civilwar 
                 + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                 ,
                 data = CR, family = binomial)

cwchal_cl2 <- coeftest(cwchal2, vcov = vcovCL, cluster = ~ cowcode)

cwchal3 <- brglm(CRbin ~ civilwar 
                 + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                 + incumbpowerdur + incgovmilitary
                 ,
                 data = CR, family = binomial)

cwchal_cl3 <- coeftest(cwchal3, vcov = vcovCL, cluster = ~ cowcode)

cwsuc1 <- brglm(sucCRbin ~  civilwar 
                ,
                data = CRsub, family = binomial)

cwsuc_cl1 <- coeftest(cwsuc1, vcov = vcovCL, cluster = ~ cowcode)

cwsuc2 <- brglm(sucCRbin ~  civilwar 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                ,
                data = CRsub, family = binomial)

cwsuc_cl2 <- coeftest(cwsuc2, vcov = vcovCL, cluster = ~ cowcode)

cwsuc3 <- brglm(sucCRbin ~  civilwar 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                + incumbpowerdur + incgovmilitary
                ,
                data = CRsub, family = binomial)

cwsuc_cl3 <- coeftest(cwsuc3, vcov = vcovCL, cluster = ~ cowcode)

## Rev militia
revmil1 <- brglm(sucCRbin ~ revcoerc 
                 ,
                 data = CR, family = binomial)

revmil_cl1 <- coeftest(revmil1, vcov = vcovCL, cluster = ~ cowcode)

revmil2 <- brglm(sucCRbin ~ revcoerc 
                 + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                 ,
                 data = CR, family = binomial)

revmil_cl2 <- coeftest(revmil2, vcov = vcovCL, cluster = ~ cowcode)

revmil3 <- brglm(sucCRbin ~ revcoerc 
                 + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                 + incumbpowerdur + incgovmilitary
                 ,
                 data = CR, family = binomial)

revmil_cl3 <- coeftest(revmil3, vcov = vcovCL, cluster = ~ cowcode)

revmilchal1 <- brglm(CRbin ~ revcoerc 
                     ,
                     data = CR, family = binomial)

revmilchal_cl1 <- coeftest(revmilchal1, vcov = vcovCL, cluster = ~ cowcode)

revmilchal2 <- brglm(CRbin ~ revcoerc 
                     + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                     ,
                     data = CR, family = binomial)

revmilchal_cl2 <- coeftest(revmilchal2, vcov = vcovCL, cluster = ~ cowcode)

revmilchal3 <- brglm(CRbin ~ revcoerc 
                     + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                     + incumbpowerdur + incgovmilitary
                     ,
                     data = CR, family = binomial)

revmilchal_cl3 <- coeftest(revmilchal3, vcov = vcovCL, cluster = ~ cowcode)

revmilsuc1 <- brglm(sucCRbin ~  revcoerc 
                    ,
                    data = CRsub, family = binomial)

revmilsuc_cl1 <- coeftest(revmilsuc1, vcov = vcovCL, cluster = ~ cowcode)

revmilsuc2 <- brglm(sucCRbin ~  revcoerc 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    ,
                    data = CRsub, family = binomial)

revmilsuc_cl2 <- coeftest(revmilsuc2, vcov = vcovCL, cluster = ~ cowcode)

revmilsuc3 <- brglm(sucCRbin ~  revcoerc 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbpowerdur + incgovmilitary
                    ,
                    data = CRsub, family = binomial)

revmilsuc_cl3 <- coeftest(revmilsuc3, vcov = vcovCL, cluster = ~ cowcode)

## Remove high-leverage observations ##############################################
plot(deaths,which = 4, id.n = 3)
# 12 "Hungarian Revolution of 1919"
# 14 "March on Rome"
# 71 "Egyptian Revolution 2011"

deaths_lev <- brglm(sucCRbin ~ lndeaths 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbpowerdur + incgovmilitary
                    + leftist + vanguard + nr_spons
                    ,
                    data = CR[-(12),], family = binomial)

deaths_lev_cl <- coeftest(deaths_lev, vcov = vcovCL, cluster = ~ cowcode)

plot(cw,which = 4, id.n = 3)
# 12 "Hungarian Revolution of 1919"
# 14 "March on Rome"
# 72 Yemeni Uprising 2011"

cw_lev <- brglm(sucCRbin ~ civilwar 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                + incumbpowerdur + incgovmilitary
                + leftist + vanguard + nr_spons
                ,
                data = CR[-(12),], family = binomial)

cw_lev_cl <- coeftest(cw_lev, vcov = vcovCL, cluster = ~ cowcode)

plot(revmil,which = 4, id.n = 3)
# 12 "Hungarian Revolution of 1919"
# 4 "Xinhai Revolution"
# 71 "Egyptian Revolution 2011"

revmil_lev <- brglm(sucCRbin ~ revcoerc 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbpowerdur + incgovmilitary
                    + leftist + vanguard + nr_spons
                    ,
                    data = CR[-(12),], family = binomial)

revmil_lev_cl <- coeftest(revmil_lev, vcov = vcovCL, cluster = ~ cowcode)

plot(deathschal,which = 4, id.n = 3)
# 98 "Ethiopian Civil War"
# 38 "Carnation Revolution" 
# 36 "Sandinista Revolution"

deathschal_lev <- brglm(CRbin ~ lndeaths 
                        + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                        + incumbpowerdur + incgovmilitary
                        + leftist + vanguard + nr_spons
                        ,
                        data = CR[-(98),], family = binomial)

deathschal_lev_cl <- coeftest(deathschal_lev, vcov = vcovCL, cluster = ~ cowcode)

plot(cwchal,which = 4, id.n = 3)
# 98 "Ethiopian Civil War"
# 38 "Carnation Revolution" 
# 36 "Sandinista Revolution"

cwchal_lev <- brglm(CRbin ~ civilwar 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbpowerdur + incgovmilitary
                    + leftist + vanguard + nr_spons
                    ,
                    data = CR[-(98),], family = binomial)

cwchal_lev_cl <- coeftest(cwchal_lev, vcov = vcovCL, cluster = ~ cowcode)

plot(revmilchal,which = 4, id.n = 3)
# 98 "Ethiopian Civil War"
# 38 "Carnation Revolution" 
# 36 "Sandinista Revolution"

revmilchal_lev <- brglm(CRbin ~ revcoerc 
                        + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                        + incumbpowerdur + incgovmilitary
                        + leftist + vanguard + nr_spons
                        ,
                        data = CR[-(98),], family = binomial)

revmilchal_lev_cl <- coeftest(revmilchal_lev, vcov = vcovCL, cluster = ~ cowcode)

plot(deathssuc,which = 4, id.n = 3)
# 12 "Hungarian Revolution of 1919"
# 100 "1990 Revolution in Bangladesh"
# 29 "Guinea-Bissau revolution"

deathssuc_lev <- brglm(sucCRbin ~  lndeaths 
                       + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                       + incumbpowerdur + incgovmilitary
                       + leftist + vanguard + nr_spons
                       
                       ,
                       data = CRsub[-(7),], family = binomial)

deathssuc_lev_cl <- coeftest(deathssuc_lev, vcov = vcovCL, cluster = ~ cowcode)

plot(cwsuc,which = 4, id.n = 3)
# 12 "Hungarian Revolution of 1919"
# 72 Yemeni Uprising 2011"
# 8 "October Revolution and Civil War (Russia)"

cwsuc_lev <- brglm(sucCRbin ~  civilwar 
                   + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                   + incumbpowerdur + incgovmilitary
                   + leftist + vanguard + nr_spons
                   ,
                   data = CRsub[-(7),], family = binomial)

cwsuc_lev_cl <- coeftest(cwsuc_lev, vcov = vcovCL, cluster = ~ cowcode)

plot(revmilsuc,which = 4, id.n = 3)
# 12 "Hungarian Revolution of 1919"
# 8 "October Revolution and Civil War (Russia)"
# 4 "Xinhai Revolution"

revmilsuc_lev <- brglm(sucCRbin ~  revcoerc 
                       + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                       + incumbpowerdur + incgovmilitary
                       + leftist + vanguard + nr_spons
                       ,
                       data = CRsub[-(7),], family = binomial)

revmilsuc_lev_cl <- coeftest(revmilsuc_lev, vcov = vcovCL, cluster = ~ cowcode)

## Alternative old regime controls #############################################################################
## Deaths
deaths_allreg <- brglm(sucCRbin ~ lndeaths 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                + incumbgovttype
                + leftist + vanguard + nr_spons
                ,
                data = CR, family = binomial)

deaths_allreg_cl <- coeftest(deaths_allreg, vcov = vcovCL, cluster = ~ cowcode)

deathschal_allreg <- brglm(CRbin ~ lndeaths 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbgovttype
                    + leftist + vanguard + nr_spons
                    ,
                    data = CR, family = binomial)

deathschal_allreg_cl <- coeftest(deathschal_allreg, vcov = vcovCL, cluster = ~ cowcode)

deathssuc_allreg <- brglm(sucCRbin ~  lndeaths 
                   + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                   + incumbgovttype
                   + leftist + vanguard + nr_spons
                   
                   ,
                   data = CRsub, family = binomial)

deathssuc_allreg_cl <- coeftest(deathssuc_allreg, vcov = vcovCL, cluster = ~ cowcode)

## Civil war
cw_allreg <- brglm(sucCRbin ~ civilwar 
            + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
            + incumbgovttype
            + leftist + vanguard + nr_spons
            ,
            data = CR, family = binomial)

cw_allreg_cl <- coeftest(cw_allreg, vcov = vcovCL, cluster = ~ cowcode)

cwchal_allreg <- brglm(CRbin ~ civilwar 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                + incumbgovttype
                + leftist + vanguard + nr_spons
                ,
                data = CR, family = binomial)

cwchal_allreg_cl <- coeftest(cwchal_allreg, vcov = vcovCL, cluster = ~ cowcode)

cwsuc_allreg <- brglm(sucCRbin ~  civilwar 
               + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
               + incumbgovttype
               + leftist + vanguard + nr_spons
               ,
               data = CRsub, family = binomial)

cwsuc_allreg_cl <- coeftest(cwsuc_allreg, vcov = vcovCL, cluster = ~ cowcode)

## Rev militia
revmil_allreg <- brglm(sucCRbin ~ revcoerc 
                + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                + incumbgovttype
                + leftist + vanguard + nr_spons
                ,
                data = CR, family = binomial)

revmil_allreg_cl <- coeftest(revmil_allreg, vcov = vcovCL, cluster = ~ cowcode)

revmilchal_allreg <- brglm(CRbin ~ revcoerc 
                    + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                    + incumbgovttype
                    + leftist + vanguard + nr_spons
                    ,
                    data = CR, family = binomial)

revmilchal_allreg_cl <- coeftest(revmilchal_allreg, vcov = vcovCL, cluster = ~ cowcode)

revmilsuc_allreg <- brglm(sucCRbin ~  revcoerc 
                   + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                   + incumbgovttype
                   + leftist + vanguard + nr_spons
                   ,
                   data = CRsub, family = binomial)

revmilsuc_allreg_cl <- coeftest(revmilsuc_allreg, vcov = vcovCL, cluster = ~ cowcode)

## Military personnel #############################################################################
## Deaths
deaths_milper <- brglm(sucCRbin ~ lndeaths 
                       + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                       + milper1000
                       + leftist + vanguard + nr_spons
                       ,
                       data = CR, family = binomial)

deaths_milper_cl <- coeftest(deaths_milper, vcov = vcovCL, cluster = ~ cowcode)

deathschal_milper <- brglm(CRbin ~ lndeaths 
                           + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                           + milper1000
                           + leftist + vanguard + nr_spons
                           ,
                           data = CR, family = binomial)

deathschal_milper_cl <- coeftest(deathschal_milper, vcov = vcovCL, cluster = ~ cowcode)

deathssuc_milper <- brglm(sucCRbin ~  lndeaths 
                          + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                          + milper1000
                          + leftist + vanguard + nr_spons
                          ,
                          data = CRsub, family = binomial)

deathssuc_milper_cl <- coeftest(deathssuc_milper, vcov = vcovCL, cluster = ~ cowcode)

## Civil war
cw_milper <- brglm(sucCRbin ~ civilwar 
                   + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                   + milper1000
                   + leftist + vanguard + nr_spons
                   ,
                   data = CR, family = binomial)

cw_milper_cl <- coeftest(cw_milper, vcov = vcovCL, cluster = ~ cowcode)

cwchal_milper <- brglm(CRbin ~ civilwar 
                       + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                       + milper1000
                       + leftist + vanguard + nr_spons
                       ,
                       data = CR, family = binomial)

cwchal_milper_cl <- coeftest(cwchal_milper, vcov = vcovCL, cluster = ~ cowcode)

cwsuc_milper <- brglm(sucCRbin ~  civilwar 
                      + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                      + milper1000
                      + leftist + vanguard + nr_spons
                      ,
                      data = CRsub, family = binomial)

cwsuc_milper_cl <- coeftest(cwsuc_milper, vcov = vcovCL, cluster = ~ cowcode)

## Rev militia
revmil_milper <- brglm(sucCRbin ~ revcoerc 
                       + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                       + milper1000
                       + leftist + vanguard + nr_spons
                       ,
                       data = CR, family = binomial)

revmil_milper_cl <- coeftest(revmil_milper, vcov = vcovCL, cluster = ~ cowcode)

revmilchal_milper <- brglm(CRbin ~ revcoerc 
                           + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                           + milper1000
                           + leftist + vanguard + nr_spons
                           ,
                           data = CR, family = binomial)

revmilchal_milper_cl <- coeftest(revmilchal_milper, vcov = vcovCL, cluster = ~ cowcode)

revmilsuc_milper <- brglm(sucCRbin ~  revcoerc 
                          + lngdpcapplus1 + lntpop + urbpercaftrev + ethfrac + lnmtnest + endyearcent
                          + milper1000
                          + leftist + vanguard + nr_spons
                          ,
                          data = CRsub, family = binomial)

revmilsuc_milper_cl <- coeftest(revmilsuc_milper, vcov = vcovCL, cluster = ~ cowcode)
